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CHAPTER I . INTRODUCTION 

Aircraft and missile fuselage geometry can often be approximated by 
circular or elliptical cones and thus designers have sought methods for 
computing the inviscid supersonic flow about these shapes. In addition, 
studies based on simple conical geometry provide a clearer insight into 
fundamental physical processes for both the experimental and computational 
investigator . 

Busemann (1) pioneered the concept of conical flow defined as a self- 
similar flov,< in which fluid properties remain constant along rays emanating 
from the conical origin. This reduces' from three to two the number of 
independent spatial dimensions in the governing nonlinear partial differ- 
ential equations. Taylor and Maccoll (2) considered the case of a 
circular cone at zero angle of attack. The flowfield, being axisymmetric, 
is determined by only one independent variable thus resulting in a system 
of ordinary differential equations. 

Inclined cones were first treated by Stone (3) as perturbations about 
the Taylor-Maccoll solutions. Using this method, which is applicable only 
for small angles of attack, Kopal (4) compiled tables of numerical results. 
Ferri (5), however, recognized that Stone ' s method was conceptually wrong 
near the surface of the cone in that it yields values for the surface 
entropy that are periodic in the meridional angle. Instead, since the cone 
surface is a streamline of the flow, the entropy must be invariant. Ferri 
amended Stone's results by showing that there exists a thin layer of 
rapidly changing entropy near the cone surface which he lobeled the vor- 
tical layer. His analysis of the crossflow velocity (the projection of 
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the velocity vector (Figure la) onto the surface of a unit sphere centered 
at the cone apex (Figure lb) ) revealed that for inviscid conical flow all 
streamlines must converge to a conical stagnation point in the leeward 
symmetry plane; the so-called vortical singularity. Since each streamline 
converging on the vortical singularity has passed through different por- 
tions of the bow shock, the entropy, density, and spherical radial velocity 
component u are multivalued at the singularity. 

Based on this theoretical background and early numerical studies a 
considerable volume of research on nonaxisymmetric conical flows has 
evolved. Numerous surveys of this work have been presented. Of these. 
Reference 6 offers the most detailed discussion of many of the methods that 
have been developed and Reference 7 contains a representative bibliography 
of the work prior to 1972. 

Until recently, however, all techniques have been limited. Even for 
the simplest conical body, the circular cone, the methods are restricted 
to cases in which the angle of attack a does not greatly exceed the cone 
half-angle 6^ or are restricted to solving only the windward portion of 
the flowfield. 

Inverse methods (8, 9) , the method of integral relations (10) , vari- 
ations of the method of lines (7, 11-14), the method of characteristics 
(13, 15) , as well as the BVLR finite-difference (16, 17) and nonconserva- 
tive finite-difference formulations (18, 19), encounter difficulties at 
large angles of attack ^ 1) • "^he nximerical problems are attributable 

to the appearance of new flow features on the leeward side of the cone. 
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Figure 
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(a) SPHERICAL POLAR COORDINATE SYSTEM AND VELOCITIES 
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VORTICAL SINGULARITY 


CONICAL STAGNATION^: 
POINTS V»W = 0 



(b) CONICAL STREAMLINES (UNIT SPHERE PROJECTION) 


1. Supersonic flow about a circular cone at large angle of attack, 

a/0 > 1 
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At moderate angles of attack the crossflow velocity is everywhere 
subsonic. At large incidence the crossflow accelerates from the windward 
symmetry plane to supersonic velocities in the shoulder region of the cone. 
In order to satisfy the leeward symmetry plane boundary conditions, an em- 
bedded conical crossflow shock (not present at smaller angles of attack) 
forms to decelerate the crossflow (Figure la, b) . 

With increasing angle of attack the bow shock gains strength in the 
windward symmetry plane while becoming weaker on the leeward side, thus 
the jump in entropy at the vortical singularity intensifies. 

The nature of the vortical singularity for cones at large angles of 
attack has been investigated in several theoretical (5, 20 - 22) and some 
recent numerical (23, 24) studies. Experimental evidence (25 - 28) indi- 
cates that an analogous feature may even be present in viscous flows. The 
theoretical predictions for the behavior of the vortical singularity are 
based upon linear theory or localized solutions and thus cannot account 
for the effects of the crossflow shock. However, the conjecture that the 
vortical singularity will, at large angles of attack, lift off the cone 
surface in the leeward symmetry plane (as shown in Figure 1) has been 
supported by the numerical results (23, 24) . 

Explicit finite-difference techniques (23, 29 - 31) utilizing the 
conservation- law form of the governing equations provided the first evi- 
dence for the inviscid, nonlinear flow pattern in the leeward region of 
highly inclined cones. The success of this approach is ascribed to the 
use of the conservation- law dependent variables. The internal crossflow 
shock is then automatically captured to within a few mesh intervals. 
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IiTtproved accuracy is obtained by treatii^g the peripheral bow shock as a 
sharp discontinuity forming the outer boundary of the computational domain, 
but the captured vortical singularity remains difficult to resolve. Fur- 
ther, the ability to capture the internal shock is bounded by an upper 
limit on the range of parameters ot/9^ and the free-stream Mach number 
beyond which the growing strength of the crossflow shock leads to 
numerical instability. This range has recently been extended by altering 
the circular shape of the cone on the leeward side (32) or through stress 
terms added to the governing equations (32, 33) . The resulting solutions, 
however, no longer represent the inviscid cone flowfield. 

An investigation preliminary to this study indicated that modifica- 
tions to the shock-capturing approach, such as the patching of the shock 
jump conditions at the crossflow shock-cone surface intersection or the 
selective use of dissipation, affords only limited improvements in resolu- 
tion and stability. 

In contrast to the finite-difference approach, a modification of the 
method of lines (24, 34) has been used to solve the leeward flow region 
provided the windward and leeward crossflow sonic lines extend to the 
bow shock (Figure 2a) . The method iterates on the internal shock location 
but utilizes a bow shock shape extrapolated from the shoulder region. The 
solutions obtained have been the most extensive to date, particularly in 
regard to the vortical singularity. The accuracy of the procedure is dif- 
ficult to assess since there is considerable disagreement with the Shock- 
captured finite-difference results. In addition, it is expected that at 
large angles of attack the possible nonmonotonic shape of the bow shock 
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Figure 2. Numerical approach determined by the extent of the zone of 

supersonic crossflow, M ^ = [ (v^ + w^)/a^]"^ 

cf 
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(35) (the maximum shock standoff distance occurring away from the leeward 
symmetry plane) and its vanishingly weak strength, will dictate a more 
accurate shock boundary treatment. 

An alternative finite-difference method which does not rely on the 
shock-capturing ability inherent in the conservation- law formulation has 
been developed in a recent series of reports (36 - 42) . Moretti refers to 
the technique as floating shock fitting to distinguish it from other sharp 
shock schemes which treat shocks as boundaries of computational domains. 

In the floating-fitting procedure shocks propagate in the computational 
domain as discontinuity surfaces. The finite-difference mesh for the 
computation of the smooth flow regions is fixed and not forced to follow 
the shocks. Topological problems associated with fitting shocks as boun- 
daries are thereby avoided. However, the differencing scheme must be mod- 
ified to prevent forming differences across discontinuities and thus 
involves unequally spaced mesh intervals. 

Questions have been raised about the stability and programming com- 
plexity of the floating- fitting method for multidimensional flows (43 - 46) 
and all of the details of the algorithm have not yet been established. 
Successful applications of the technique, however, have demonstrated the 
feasibility of the floating- fitting approach. In particular, the method 
has been tested on the cone at zero incidence problem (37) and gave 
excellent results. 

This present study describes a numerical method which applies the 
concept of floating shock fitting to the cone at large incidence problem. 
The bow shock, the embedded crossflow shock, and the vortical singularity 
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are treated as sharp discontinuities which float simultaneously through 
the computational mesh. Use is made of the conservation-law form of the 
governing equations in the flowfield interior to aid in the detection and 
monitoring of the crossflow shock. The fitting of the crossflow shock 
avoids the stability problems encountered with the shock-capturing ap- 
proach. Further, in contrast to the numerical method of References 24 and 
34, shock layers with limited zones of supersonic crossflow adjacent to 
the cone surface (Figure 2b) can be computed. The present method is 
formulated to be applicable over the complete range of parameters M , a, 

CO 

and 0 for which the flow remains conical (bow shock wave attached) . 

c 

Along with the floating-fitting method, a technique in which the bow 
shock forms a boundary of the computational mesh has been developed. Aside 
from serving as a check on the floating-fitting approach, this "shock-as- 
a-boundary" code provides a convenient means of supplying initial condi- 
tions for the floating algorithm. 

Several boundary condition schemes have been tested in conjunction 
with the floating-fitting method. Comparisons of these floating-fitting 
solutions with shock-as-a-boundary results are presented together with 
numerical results obtained in References 13, 31, and 34 and with experi- 
mental measurements (47) . 
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CHAPTER II. GOVERNING EQUATIONS — CHOICE OF VARIABLES 

The selection of appropriate dependent and independent variables can 
simplify the application of a numerical technique and result in improved 
accuracy. For example, enhanced shock-capturing properties are obtained 
through the choice of the conservation-law form of the dependent variables 
in a coordinate system that closely parallels the shock. The implementa- 
tion of boundary condition schemes, however, may often best be achieved 
through the use of the primitive dependent variables in a boundary oriented 
coordinate frame. 

Several forms of the governing equations, each chosen with regard to 
its application and the flow region being computed, are used in the cone 
at angle of attack problem. 

The equations governing the flow of an inviscid, nonconducting, com- 
pressible fluid in spherical polar coordinates (Figure la) may be arranged 
into weak conservation- law form as 

+ E^ + F- + G , + H = 0 (1) 

t R D 9 

where U, E, F, G, and H are the five- component vectors, 

pu 

p + pu^ ^ 

E = ■* p uv f , F 

puw 

(e+p)u 


p (2u + V cot 0) 
p ( 2u^ - w^ - v^ + uv cot 0 ) 
H = — p [3uv + (v^ - w^) cot 0] 

^ pw(3w + 2v cot 0) 

(e + p) (2u + V cot 0) 
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For an ideal gas the system is completed using the equation of state to 
relate the total energy per unit volume e to the pressure p, density p, 
and velocity components u, v, and w by 

e = p/(y- 1) + p (u^ + v^ + w^) /2 (3) 

where y ratio of specific heats. The pressure and density are 

referenced to free-stream static conditions. The velocity components are 
thus normalized by where a^ is the free-stream speed of so\ind. 

If the flow is assiamed to be steady (U. = 0) , the integrated form of 

t 

the energy equation replaces the last row of Equation (1) with the alge- 
braic relation 

p = p{ [1 + (Y - l) M^/2] - [ (y - 1) + w^) ] /2 y) (4) 

If the flowfield is further assumed to be conical (E = 0) , the govern- 
ing equations in the self-similar crossflow plane (R=l) become 

Fj + + H - 0 (5) 

Equation (5) is elliptic for cones at small or moderate angles of 
attack. With increasing incidence a small zone of supersonic crossflow 
forms adjacent to the surface on the side of the cone (Figure 2b) . At 
still larger angles of attack, the region of supersonic crossflow extends 
to the bow shocks This enables partitioning of the crossflow plane into 
Windward and leeward elliptic regions separated by the hyperbolic shoulder 
region (Figure 2a) . 

In the application of the floating-fitting method. Equation (5) is 
rendered hyperbolic by the addition of the unsteady derivative term. 
Alternatively, in the algorithm which fits the bow shock as a boundary, 
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the spherical radially-asymptotic approach (addition of the conical 
derivative term) is followed. The unsteady approach is pi'eferred since 
the hyperbolic nature is preserved over the entire range of parameters 
M , a, and 9 for which the bow shock remains attached. Thus, unlike 

CO c 

the radially-asymptotic method, the unsteady approach allows computation 
of radially subsonic flows (u < a, where a = local speed of sound) which 
can occur at large angles of attack (23) . 

An a priori assumption (based upon previously obtained solutions) can 
often be made regarding the size of the region of supersonic crossflow. 
When partitioning of the crossflow plane is possible (Figure 2a) , maximum 
efficiency is achieved by solving each region separately. In addition, 
for a given number of computational plane grid points resolution is 
improved since the mesh point aensity in each region increases. More 
importantly, the leeward region is freed from the small step size imposed 
by the stability condition (presented in Chapter III) in the windward 
symmetry plane. Another benefit of partitioning is that a substantial 
part of the shoulder region can be computed in one less dimension by a 
(j>-direction marching code which integrates the crossflow plane equations. 
The shoulder region algorithm (described in Appendix A) starts from the 
windward region outflow boundary (with w > a) and sweeps around the cone 
until the meridional Mach number approaches unity. 

Reference MeSh 

The floating-fitting procedure utilizes a fixed reference mesh de- 
fined by (j) = constant meridional boundaries and a circular outer boundary 
An independent variable transformation 


9 = 9 
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X 


0 - e ^ 

c 


9-0 
o c 


Y = } 


( 6 ) 


T 


= t 


normalizes the cone-to-outer boundary distance. The reference mesh for 
the leeward region calculation is illustrated in Figure 3. 

The transformation Equation (6) distinguishes the floating-fitting 
technique from shock-as-a-boundary methods (cf . Equation (A4) ) and reduces 
the complexity of the flowfield interior tr,..n&formed governing equations. 
This simplification is important in the leeward region where the shock-as- 
a-boundary approach requires normalizations in two directions. 

A further transformation (X, Y, x) , which permits clustering points 
near the cone surface and about a chosen meridional location (48) , may be 
used to concentrate points near the crossflow shock. Clustering is also 
useful in keeping the number of reference mesh points that lie outside of 
the shock layer (where no calculations take place) to a minimum. This is 
especially helpful when solving the symmetrical half-plane (Figure 2b) 
problem since the bow shock is often far from being circular. The clus- 
tered reference mesh independent variables are defined by 


X = X(X) 
Y = Y(Y) 
T = f ( t) 


= — sinh" ^ [X sinh 


max 


B + sinh" ^ 


sinh B 


(7) 


where 








14 


$ = Y - 6 

Jj 

$ = tf) - 6 

max u L 


B = -ln 


1+ (e’^- 1) 

$ 

o 

^max 


1 + (e"'^ - 1) 


^max^ 



( 8 ) 


with 


(j) = lower meridional boundary 

L 

(f) = upper meridional boundary 

ij; = meridional clustering parameter, clustering about (j)^ 
ui = theta clustering parameter, clustering about 6 


Application of the clustering transformation (Equation (7)) to the 

conical Euler equations in spherical polar coordinates (Equation (1) with 

S =0) and rearrangement into weak conservation-law form yields the system 
R 


U_ + F- + G- + H = 0 
T X Y 


(9) 


where 


U = U 


— 1 — rv/ 

o c 


( 10 ) 


G = y^G 


H = H 


(Xy)^F - 


6 -6 
o c 
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The first derivatives of the clustered mesh independent variables appearing 
in Equation (10) are 


= 


sinh w 


X (jj cosh(wX) 


y.. = 


1 


u 


sinh B 


r ■ 

f- j 

1 

sinh B 

, 

+ 1 

. 

1 $ j 

J 



( 11 ) 


and the second derivatives required to form the H vector are 


Wx)x 


(Vy)- - 


sinh (u sinh(g)X) 
cosh^ (wX) 


$ 


- 1 


sinh^ B 


ii + 

$ 

1 

sinh B 

L 

$ 

o 



( 12 ) 


The number of reference mesh points in the X and Y--directions, the 
outer boundary location 6^, and the clustering parameters ijj , and cj>^, 

are specified in a separate code as part of the initialization procedure 
for the floating- fitting code (initial conditions are discussed in Chapter 
IV) . This mesh generation code offers an interactive option in which the 
mesh variables may be iteratively adjusted. The resulting distribution of 
grid points in the physical domain is displayed on a cathode ray tube 
(C.R.T.) following each mesh specification (IBM 360/67 time-sharing system 
linked to an IBM 2250 C.R.T. ). 


Once an initial reference mesh is selected, it remains fixed in the 
unsteady relaxation procedure. Thus the geometric derivatives, which 
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appear as coefficients of the F and G vectors in Equation (10 ) , need 
only be evaluated initially for each mesh point. 

Floating Discontinuity Points 

The floating- fitting technique also employs discontinuity oriented 
coordinate systems, each normalizing the distance between a fixed boundary 
and a moving discontinuity (Figure 3) . The necessary transformations are 
presented in Table 1. 

The governing equations are integrated on one side of each discontin- 
uity. At embedded shock discontinuities, the integration is performed on 
the low pressure side. At the vortical singularity, the integration is 
performed on the low entropy side (either side may be used) . The low 
pressure side integration, however, is not required at the peripheral bow 
shock since the flowfield is the known free stream. The flow values on 
the remaining side of each discontinuity are obtained, once ti.e new dis- 
continuity geometry has been determined, from the discontinuity jump con- 
ditions. 

The form of the Euler equations suitable for all discontinuity 
alignment transformations on the unit sphere is 

d- + [Bi] d- + [Bz] d- + A3 = 0 (13) 

where 


r > 

p 


a^p (2u + V cot S) 

u 


-(v^+w^) 

V 

A3 = ■ 

uv - w^ cot 6 

w 


w (u + V cot 0 ) 

.e. 


_ (e + p) (2u + V cot 9) 


d = ( 


(14) 
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Table 1. Discontinuity alignment transformations 


Coordinate system 


6-direction 


c() -direction 


0 - 

9 


c 

0 - 

6 

0 

c 

0 - 

6 


c 

9 - 

9 

s 

c 

X 


X 


0 - 

0 


0 

9 - 

9 

0 

vs 

0 - 

9 


c 

9 

- 6 

vs 

c 


Reference mesh 


Bow shock, high-pressure side 


a 


Crossflow shock, low-pressure side' 




Crossflow shock, high-pressure side 




Vortical singularity, bow shock side 


Vortical singularity, body side 


d 


X = 


Xl = 


X2 = X 


X3 = X 


X4 = 


X5 = 


Y 


Yi = Y 


^2 = 


Yq = 




(j) - (j> 


^ <)> - 

u ^s 


Y4 = Y 


Y5 = Y 


a 


9 =6 (<J),t) ~ bow shock shape 

s s 


(j)^ = (()^(6,t) ~ crossflow shock shape 
s s 


~ upper meridional boundaries 


^vs ~ ^vs^^^ ~ vortical singularity position 
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with 


[Bi] = {Xq[Ai] + X [A2] + X^[I]}X, 


X 


[B 2 ] = {Yg[Ail + Y^[A2] + 


( 15 ) 


where [I] is the identity matrix and 


[All = 


[A 2 ] = 


V 

0 a 

2p 

0 

o' 



0 

V 

0 

0 

0 



1/p 

0 

V 

0 

0 

/ 


0 

0 

0 

V 

0 



V 

0 e 

+ p 

0 

V 





w 

0 

0 

a 

2p 

0 ’ 

1 


0 

w 

0 


0 

0 

-L 


0 

0 

w 


0 

0 

sin 6 



1/P 






0 

0 


w 

0 



w 

0 

0 

e + p 

w 


(16) 


The derivatives of the reference mesh independent variables are evaluated 
for each alignment transformation i=l, 2, 3, 4, 5 (Table 1) by applica- 
tion of the chain rule where 

- !<e0Xi% 

= X0 0x. Xi^. 

^0 " ^4> ^i0 

\ ’'it 


(17) 


with Xq = 1/(6 - 0 ) and Y = 1. 
y o c (p 

To integrate Equation (13) , the reference mesh derivatives (listed in 

Table 2) are evaluated at each discontinuity. The coefficients of: 0o./ 

<P 

6g^, 0vst' *^30' Table 2, then reduce to -1^ 

The nonconservative form of the governing equations (Equation (13)) 


does not require evaluation of the second derivatives of the transformed 




Table 

2. 

Reference mesh 

independent variable derivatives 




Coordinate System 



X-direction 




Y-direction 




^9 





^9 



^t 


Spherical polar 

1 


0 

0 


0 


1 

0 


Reference mesh 

1 


0 

0 


0 


1 

0 


9 - 
o 

9 

c 




Bow shock , high- 

1 


0 

-X 


0 


1 

0 


pressure side 

9 - 
o 

9 

c 

9-9 

sc 

9-9 
s c 




Crossflow shock, 

1 


0 

0 


-(V-y 




d)~ 

low-pressure side 

9 - 
o 

9~ 

c 


tl 

-e- 

1 

-e- 



Crossflow shock. 

1 


0 

0 


-(Y - (j) ) 
u 



■«- V 


high-pressure side 

9 - 

o 

9 

c 


^s u 

J. 

(l)~- (f) 
s ^u 

Vortical singularity. 

1 


0 

-(X- 1) 

®vs^ 

0 


1 

0 


bow shock side 

9 - 

o 

9 

c 

9 - 9 

vs o 




Vortical singularity. 

1 


0 

-X 

9vst 






body side 

9 - 
o 

9" 

c 

9 - 9 

vs c 

0 


1 

0 
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independent variables. These terms, which contribute to the source vector 
of the conservation-law system (see, for example. Equation (10)), are 
absent in the nonconservative system. 

The choice of circular geometry for the reference mesh is responsible 
for the simple scaling which relates derivatives along discontinuities and 
derivatives in the discontinuity normalized direction to the reference mesh 
system. The finite-difference representation of Equation (13) for a dis- 
continuity integration does not require grid points located at AX^ and AY^ 
(Table 1) intervals. Each discontinuity mesh consists solely of points 
lying on the discontinuity and being tracked along one of the reference 
mesh directions (Figure 3) . 
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CHAPTER III. FINITE-DIFFERENCE METHOD 

MacCormack ' s explicit, second-order, noncentered, predictor-corrector 
algorithm (49) (presented in Appendix B) is used to advance the flowfield 
in the time- asymptotic relaxation process. The algorithm's low storage 
requirements and ease of programming make it readily adaptable to a 
floating-fitting approach. With floating-fitting, the standard MacCormack 
scheme is modified for mesh points which lie on the boundaries and float- 
ing discontinuities; and for points neighboring (in time or space) the 
floating discontinuities. In the reference mesh interior the forward- 
predictor, backward-corrector sequence is followed provided the computa- 
tional molecule is not crossed by a discontinuity. 

Special Discretization Formulas 

The floating-fitting technique introduces Unequally spaced mesh inter- 
vals. Special differencing approximations (36) , formulated to maintain 
stability and accuracy, replace the two-point uniformly spaced approxima- 
tions used at each level in the MacCormack algorithm. The derivation of 
these differencing approximations is based upon weighting the contributions 
of neighboring mesh points (on the same side of a discontinuity) in such a 
way that the truncation error varies smoothly as the discontinuity cuts 
through the mesh. 

Reference mesh points 

At a given time level, the presence of a discontinuity adjacent to a 
reference mesh point will prevent forming the standard spatial difference. 
In addition, the temporal derivative must be discarded if, in advancing 
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to the next level, the reference mesh point passes from one side of the 
discontinuity to the other. 

Temporal derivatives The floating-fitting algorithm keeps track 

of the location of moving discontinuities relative to the fixed reference 
mesh. Movement of discontinuities over reference mesh points generally 
occurs only during the early transient phase. Consider the integration of 

U +F +G„+H=0 (18) 

T X Y 


where 



(19) 


H = H 


using MacCormack's scheme. Referring to Figure 4, if the shock point S 
is located between reference mesh points J and J-1 at time level n and 
crosses point J (moving in the X-direction) while advancing to time level 

H*f'l I 

n+1, then U is computed by the linear interpolation 
J 


.n+1 


1 + e 


X 




( 20 ) 


Spatial derivatives The replacement of reference mesh spatial 
derivatives is illustrated in Figure 4. The approximation (Salas, M. , 
private communication, NASA- Langley Research Center, 1976) shown for F^j 
was found to yield results equivalent to those obtained using a formula 
that involved an additional point in the reference mesh (37) . The removal 
of even one point offei^s coding simplifications. 



FLOATING 
SHOCK POINTS 

(X-DIRECTION INTERPOLATED 
TRACKING) SHOCK POINTS 



REFERENCE MESH 
COMPUTATIONAL 
PLANE 


DIFFERENCE APPROXIMATION FOR ~ 
AT SHOCK POINT S: 


AT REFERENCE MESH POINT J : 

IN THE FORWARD PREDICTOR STEP. 
REPLACE 

*^J+l~Fj 
- -AX~ 

BY 

I r2(2-€x) 


Figure 4. Spatial differencing approximations for floating- fitting 
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If the shock in Figure 4 also cuts the mesh interval J to J-1, so 
that Fj ^ is not available for the computation of then is 

approximated by 


F - F 
S SL 


X.T - 


X„ - X 


SL 


( 21 ) 


where S-and SL represent shock points lying above and below the reference 
mesh point J. 


If point 


then is 


J is on an X-boundary, 
computed using 


so that no points are below it, 


•X- 




for 


> i 


( 22 ) 


for < 1/2, Fj^j is obtained by extrapolation of F^, in the Y- 

direction. 


Similar discretization formulas apply to the required Y-derivatives 
at point J. However, since shock points are tracked only along Y = con- 
stant lines (for the example shown in Figure 4) interpolation is required 
to obtain G_ values where the shock cuts a Ay mesh interval. Linear 
interpolation has been found to be sufficient. 

If the independent variable clustering transformation (X, Y, t) (Equa- 
tion (7)) has been applied to the reference mesh, then interpolations are 
performed in the Unequally spaced physical plane. The conservation- law 
vector in the computational plane Gg (Equation (10)) is obtained by 
multiplying the interpolated G vector by the clustered mesh Y geo- 
metric derivative. The geometric derivative is coniputed analytically 
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(Equation (11) ) using the interpolated location of the intersection of the 
shock with the physical plane A(|) mesh interval. Interpolated distances 
are transformed to the computational plane and the appropriate special 
differencing formula applied. 

Discontinuity mesh points Derivatives along the shock in Figure 4 
are computed using the usual MacCormack scheme. Derivatives in the shock 
normalized direction, obtainable in terms of X-derivatives , are computed 
using the Equation (42) shown in Figure 4. 

The same derivative approximations are used at the crossflow shock and 
at the bow shock except that the roles of X and Y are interchanged. Bow- 
shock points are tracked in the X-direction while crossflow shock points 
move in the Y-direction. 

X-derivatives at the vortical singularity are similar to those at bow 
shock points since the vortical singularity is also tracked in the X- 
direction. For the Y-derivative at the vortical singularity (w^ being 
the only nonzero Y-derivative) , interpolated reference me 'h values must 
be used. 

Stability Analysis 

Although a two-dimensional analysis is appropriate, one-dimensional 
amplification matrix theory yields a larger, yet stable, estimate for the 
step size. The influence of the discontinuity alignment transfor-mations is 
reflected by the derivatives of X and Y (Equation (17)) appearing in the 
C-F-L condition (Courar.c-Friedrichs-Lewy) 


At = min (At-, At-) 

A X 


( 23 ) 



where 


^ (CN)AX 


Ax- 


max 
(CN) AY 


^ l%ax‘'®2=> 


( 24 ) 


with 


([Bi])l 
max ^ ' 

•= max|a( [BiD 1^ 

([B2l)| 
max ' 1 

= max|a ( [B 2 ] ) | ^ 


(25) 


where i = 1, 2, 3, A, 5 (Table 1), and with the projections of the char- 
acteristic slopes (eigenvalues) determined by 


cf([Bi]) = X 


a([B2l) i = Y^ 


X^ w 

I X V + — 

6 srn 9 


Y w 

I Y„ V H r—r + Y I 

6 sin 0 t ' 


y2 4. 

X 


> 


_sin 9_, 

- 


y2 1 

f 1 

<|) 



^9 + 

sin 0_ 




(26) 


The Courant parameter, CK in Equation (24), by linear analysis must be 
1. For most floating- fitting applications CM is set to 0.9. 


Accuracy and Convergence 

In the unsteady analysis, the total enthalpy (which should approach 
that of the free stream in the converged solution) serves as a convenient 
measure of accuracy. For example, for a typical leeward region calcula- 
tion with a reference mesh containing 17 points in the X-direction and 19 
points in Y-direction the total enthalpy at interior mesh points differs 
from the free-stream total enthalpy by much less than 0.1%. Larger errors 
generally are found on the cone surface or near the vortical singularity 


(results are presented in Chapter VI) , However, the total enthalpy error 
at these boundary points stays well below 1%. 

The total enthalpy is also used to verify the achievement of a steady 
state by testing for negligible change between current values and those 
computed 10 iterations before. Due to the extensive use of interactive 
graphics and a time-sharing system, meaningful convergence times are dif- 
ficult to obtain. Roughly, a 17 x 19 mesh requires less than 15 minutes on 
an IBM 360/67 to achieve convergence. 
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CHAPTER IV. INITIAL CONDITIONS 

Free-stream conditions or very approximate initial guesses are often 
sufficient to initiate a shock-capturing relaxation procedure. Somewhat 
more care is required in specifying the initial flowfield when a peripheral 
shock is computed by a shock fitting approach. If all discontinuities are 
to be fitted, the determination of initial conditions can become very in- 
volved. In the case of a cone at large angle of attack the problem is 
further complicated by the difficulty of obtaining a physically relevant 
flowfield that can evolve into the desired conical flow. For many cases 
of interest, the idea of starting with a sphere-cone solution is ruled out 
by the failure of existing sphere-cone codes at large angles of attack. At 
smaller angles of attack, where solutions can be obtained, the presence of 
strong entropy gradients from the blunt nose may require special treatment. 

The initialization procedure for the floating- fitting algorithm is 
illustrated schematically in Figure 5. Essentially, the strategy is to 
start from a flowfield obtained by choosing shock shape and surface dis- 
tribution parameters or by using previously computed solutions with nearly 
the same values of 9^, and a. This starting flow is then further 

refined by iterating with the code which fits the bow shock as a boundary 
and captures embedded discontinuities. Only a few iterations are required 
to establish a flowfield that is suitable for the floating-fitting tech- 
nique to adequately detect the forming embedded discontinuities. Insta- 
bilities that may arise with the capturing code, if strong crossflow shocks 
are present, are avoided by beginning the floating-fitting technique in tlie 
early stages of shock formation. 
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Figure 5. Determination of initial conditions for the floating- fitting scheme 
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Bow Shock Shape 


As shown in Figure 5, initialization of the windward and leeward re- 
gions requires an estimate for the bow shock shape. Some empirical cor- 
relation formulas have been proposed for the windward region bow shock 
shape (see, for example, Reference 9), however, a more convenient estimate 
is provided by a scheme suggested in Reference 18, The shape of the bow 
shock, of which only the windward portion is used, is assumed to be an 
off-centered ellipse when projected onto the cylindrical (r, (j)) pilane per- 
pendicular to the cone centerline (the z-axis) . In the cylindrical polar 
coordinate system, as adopted in Appendix A, tj) is measured from the wind- 
ward symmetry plane. Denoting the bow shock cylindrical radius as r^ 
and the cone surface as r^, the initial bow shock shape (at z = l) is 


r^(*) = 


/ (B^ - C^) 9 9 

-C cos (j) + By — — -9 sin“^ (p + cos'^ (p 




A 2 

sin^ (f) + cos^ (j) 


(27) 


where the parameter C is the distance in the symmetry plane that the 
center of the ellipse is offset from the center of the coordinate system 
(r = 0) . The parameters A and B are, respectively, the semiminor and 
semimajor axes of the ellipse. In terms of the dimensionless shock stand- 
off distance in the windward symmetry plane 6 , at the cone shoulder «S_„, 
and in the leeward symmetry plane 6 ^, the ellipse paraimeters are 
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B = r 


1 + 


6^+5 
L u 


C = 


A = 


~(<S -6 ) 

2 u L 

G SH 


1 - 


T 


( 28 ) 


where the standoff distances are normalized by the cone radius. 


The number of parameters that must be guessed is reduced from three to 
two by the introduction of the empirical correlation (18) , known to be 
valid over a wide range of Mach numbers and angles of attack (though un- 
tested at angles of attack exceeding the cone half-angle) , 


tan”^[(l + 6 ) tan 9 ] = 26 
u c s 


- tan"^ [ (1 + 6 ) tan 8 ] 

SH c 


a=0 


( 29 ) 


where 0 is the bow shock angle corresponding to the zero angle of attack 
s 

Taylor-Maccoll (2) solution. 

The leeward region initial bow shock may be determined without guess- 
ing shape parameters. A quadratic extension of the shoulder region bow 
shock is completely determined by the conditions of matching the known 
shock slope and standoff distance at the lower meridional boundary and by 
satisfying the zero slope condition at the leeward symmetry plane. For 
some cases, this quadratic shock may intersect the free-stream Mach cone. 

If a test reveals that bow shock mesh points lie vrithin the free-stream 
Mach cone, an incremental shock standoff distance is added. The bow shock 
is moved outward a trial amount in the leeward symmetry plane- The incre- 
mental function is specified as a cubic in <}i satisfying zero slope con- 
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straints at both leeward region meridional boundaries. The function de- 
creases from the trial value in the leeward symmetry plane to zero at the 
lower meridional (shoulder) boundary. 

Distribution of Plow Variables 

As indicated, the difficult task of specifying an initial flowfield 
for the floating- fitting code that may include embedded discontinuities is 
alleviated by first iterating on a rough initial flowfield approximation 
utilizing the capturing capabilities of the shock-as^a-bovmdary code. 

The procedure used to obtain the rough initial flowfield approximation 
may be illustrated by considering the windward region problem. Here the 
nondimensionalizations are those used in the shock-as-a-boundary code and 
described in Appendix A. 

At the bow shock, the analytically prescribed bow shock shape together 
with the known free-stream conditions and an assiamed zero shock speed are 
sufficient to determine all downstream shock values. The appropriate shock 
jump conditions are given in the Shock Boundary section of Appendix A. 

At the cone surface, a guess of the windward stagnation point pressure 
p^(0) is the only parameter that is required. Based on previous windward 
solutions, the pressure decrease around the cone may be approximated by 

where p^ is the free-stream pressure. 

The surface entropy is constant and assumed to equal that behind the 
bow shock in the windward symmetry plane, denoted as S^. Thus, surface 


(P^(0) +P„ 


) + 


(p^(0) -p^)cos 


(30) 
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density is given by 


P,(^) = 


Pc(<^) 


(31) 


The velocity modulus is obtained from the condition of constant total 
enthalpy as 




(32) 


The windward region meridional outflow boundary (|)^ is assumed to be 
positioned downstream of the crossflow sonic line. On the cone surface the 
w-velocity component must satisfy the supersonic outflow condition as well 
as the symmetry plane condition that w = 0 at = 0. 

A linear increase in w is assimed according to 


w ((jj ) = 
c 


(Y - 1) P(^*) ^ 
2 p (<}).) (j) 


(33) 


where the crossflow sonic line intersects the body at 4> = (j)^ which is 
specified as ij)^ = 

Combining constant total enthalpy (Equation (32)) with the tangency 
condition yields the remaining cylindrical polar velocity components (Ap- 
pendix A) 


Uc(4.) 




w (<t>) 
c ^ 


1 + r2 


(34) 


and 
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v^(<P) = u^((j)) rc^ (35) 

Linear interpolation between the surface and bow shock furnishes the 
interior flowfield pressure, density, w- and v-velocity components with the 
u- velocity component calculated to satisfy total enthalpy. 

In the windward symmetry plane the linear density distribution is re- 
placed in order to maintain the streamline entropy S^. The symmetry plane 
linear pressure distribution is accepted with the density determined as in 
Equation (31) . 

A similar procedure is followed to obtain a rough initial flowfield 
approximation for the leeward region problem (Figure 5} . The main differ- 
ence with the windward region technique is that the inflow meridional plane 
flow values are known. The leeward interior flowfield is computed as an 
average of a linear interpolation between assumed shock and surface dis- 
tributions and a linear interpolation in the meridional direction between 
assumed symmetry plane values and known inflow conditions. 
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CHAPTER V. BOUNDARY CONDITIONS 

In the shock-as-a-boundary approach, which captures embedded discon- 
tinuities, the physical plane bow shock layer boundaries are mapped onto 
rectangular computational plane boundaries. The flow tangency condition is 
imposed on the impermeable cone surface and symmetry planes. As a part of 
the relaxation process, use is made of the asymptotic condition that the 
surface and symmetry plane entropy can vary only through discontinuous 
jumps at the vortical singularity and at the crossflow shock-cone surface 
intersection. Of the permeable boundaries, the shock-as-a-boundary bow 
shock is governing by the Rankine-Hugoniot shock jump relations while the 
meridional boundaries encountered in the windward, shoulder, and leeward 
region subproblems are either known inflow boundaries or boundaries where 
the flow is assumed to be outflowing supersonically. 

The floating- fitting approach is distinguir,.hed from the capturing ap- 
proach by imposing jump conditions across discontinuities that do not form 
reference mesh computational plane boundaries (Figure 3) . The floating 
vortical singularity is an impermeable point in the crossflow plane projec- 
tion while the floating bow and crossflow shocks are permeable surfaces 
with, respect ' vely, known and unknown upstream flow conditions. 

One-Sided Differencing 

The implementation of boundary condition schemes at reference mesh 
points on the cone surface and the symmetry plane involves the evaluation 
of derivatives normal to these surfaces. Similarly, at floating mesh 
points on each side of a discontinuity, derivatives in the discontinuity 
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normalized direction are required. The normalized direction, as defined in 
Table 1, is not necessarily the normal direction since the alignment trans- 
formations are nonorthogonal . 

At the reference mesh boundaries the normal direction derivatives are 
evaluated using equally spaced, three-point, one-sided approximations in 
both the predictor and corrector steps of MacCormack's algorithm provided 
a discontinuity is not encountered. The use of three-point difference ap- 
proximations is in keeping with the second-order of the algorithm. Suit- 
able one-sided approximations at reference mesh boundary points in the 
vicinity of discontinuities, and for floating mesh points, have been given 
in the Special Discretization Formulas section of Chapter III. 

Cone Surface 

The shock-as-a-boundary code employs Abbett's (50) Euler predictor/ 
simple wave corrector procedure at the cone surface. The application of 
the scheme for the shoulder region problem is presented in Appendix A. 
Reference 31 outlines a version of Abbett's scheme similar to that adapted 
in the shock-as-a-boundary, radially-asymptotic, windward region code where 
the surface entropy is assigned the value behind the bow shock in the sym- 
metry plane. 

A modification to Abbett' s scheme has been devised to account for the 
jump in entropy at the crossflow shock occurring in the leeward and sym- 
metrical half-plane problems. The technique falls short of the accuracy 
obtainable with a properly formulated shock fitting procedure but is a 
vast improvement over the isentropic Abbett scheme. The modified scheme 
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also eliminates the entropy oscillations neighboring the crossflow shock 
that occur using the nonisentropic Abbett scheme presented in Reference 51. 

In the radially-asymptotic, shock-as-a-boundary approach, the modifi- 
cation to Abbett's scheme is applied following the corrector step of the 
nonisentropic Abbett procedure (51), Denote these corrected surface flow 
values as p, p, u, v, and w. Here the velocity components are in the 
cylindrical polar coordinate (z, r, tj>) directions and the nondimensionaliza- 
tions are as presented in Appendix A. Let the modified surface flow values 
be denoted by p, p, u, v, and w and let the subscripts 1 and 2 designate, 
respectively, the crossflow shock upstream and downstream conditions. 

The pressure from the Abbett scheme expansion or compression is not 
altered by the modification procedure, that is, at the z step, 

P((j)) = P(<l>) (36) 

This surface pressure distribution is scanned atd the crossflow shock is 
assumed to lie within the mesh interval containing the maximum compressive 
pressure gradient. 

Least-squares fit polynomials (2nd-degree) are passed through the 
pressure data lying on the crossflow upstream and downstream sides of this 
mesh interval providing approximations for the pressures at the crossflow 
shock of the form 

Pi (())) = aij)^ + b(p + c (37) 

P 2 (<j>) = d^^ + e<j) + f (38) 

The crossflow upstream density distribution follows, once the upstream 
side has been identified, by imposing on the cone surface the value of 
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entropy obtained at the bow shock in the windward symmetry plane, denoted 
as Sj. Thus, 


P (<j>) 


P(<|)) 

I Si J 


(39) 


The velocity components, which satisfy surface tangency after the turn 
in Abbett's scheme, are proportioned to satisfy constant total enthalpy by 
the relations 


u = u 


V = V 


where the velocity moduli are 
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(40) 


q = y 1 - 
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(41) 


P J 


The upstream surface crossflow Mach number 


M 


w 


cf 


2 p 


(42) 


is computed at each mesh point and a least-squares fit quadratic approxima- 
tion is obtained for the upstream Mach number at the crossflow shock as 

Mcf ((!>) = q<{>^ + h(|) -f i (43) 
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The location of the crossflow shock within the mesh interval contain- 
ing the maximiom compressive pressure gradient is computed by solving for 

that (() = <f>^ which satisfies the normal shock Rankine-Hugoniot jump condi- 
s 


tion 


P2(<|») - Pi (({)) 


(.j)) - (y -1) 


7 + 1 


(44) 


Equation (44) is solved iteratively using the Newton- Raphson technique 
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where m is the iteration index and 


= (y + 1) (dtj)^ + e(j) + f ) + (a^^ + b(j) + c) [ (y - 1) - 2y (gc|)^ + h<}> + i) (46) 

The mesh interval mid-point serves as the initial guess for 

With 4)^ determined, the crossflow shock entropy jump is given by 
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The crossflow downstream modified surface density and velocity distri- 
butions are computed, in the same manner as on the upstream side, but with 
the value of entropy Sj being imposed. 

The floating-fitting code uses Kentzer ' s surface boundary conditions 
(52) in conjunction with the MacCormack algorithm in a two-step procedure 
analogous to that proposed in Reference 19. Also, included as an option 


40 


in the floating- fitting code, is a simple boundary condition scheme termed 
a one-sided differencing approach. 

In the one-sided differencing scheme the conservation- law form of the 
gas dynamic equations (Equation (9)), with the tangency condition v = 0, 
is integrated on the cone surface. The resulting values for total energy 
e and the velocity components u and w are accepted. Recall that in the 
floating-fitting, time-asymptotic code the velocity components u, v, and w 
are with respect to a spherical polar (R, 6, (}>) coordinate system and the 
nondimensionalizations are as described in Chapter II. Surface pressure is 
computed by iteratively solving (using Newton-Raphson) the nonlinear equa- 
tion relating p to the total energy e 


where 


e = bp + cp 
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(48) 
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The value of entropy S appearing in the c coefficient is set to the 
appropriate crossflow shock upstream or downstream value. This surface 
entropy at the new time level is determined in a preliminary calculation 
which advances the floating bow shock only in the windward symmetry plane 
(not required if solving just the leeward region) and the floating cross- 
flow shock only at the cone surface. With p and S known, a new density 
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is obtained (see Equation (39)) and used to recompute all but the last row 
of the U vector (Equation (2) ) . 

Kentzer's surface boundary scheme is an approximation of exact charac- 
teristics theory which eliminates the interpolations and iterations 
required in following characteristic directions in a fixed mesh. 


The- approach combines the differentiated surface tangency condition 
(v_ = 0) with the compatibility relation along the down-running character- 
istic in the X, x-plane resulting in a differential equation for the sur- 
face pressure 
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With V = Of the remaining velocity components u and w are computed using 
respectively the R and (j>-momentum equations, the second and fourth rows 
of Equation (13) 
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Derivatives along the cone surface (Y-direction) are approximated by 
the MacCormack forward-backward sequence or, if neighboring a discontin- 
uity, by special discretization formulas (Chapter III). The only exception 
being the w- derivative in Equation (52) which is always approximated by 
backward differences (19) . 
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As in the one-sided scheme, the surface entropy at the new time level 
is obtained in a preparatory step by advancing the floating shocks. The 
updated density is computed with the specified entropy and the pressure 
from Equation (39) . Equation (13) , for the total energy, need not be inte- 
grated, instead, e is obtained from Equation (3) . 

Kentzer's scheme is modified for surface mesh points in the symmetry 
planes (which are conical stagnation (v = w=0) points). The symmetry con- 
ditions 



are applied and thus Equation (52) is discarded. Note that u cannot be 
obtained directly by integration since, upon substitution of w = 0, Equa- 
tion (51) reduces to 

u_ = 0 (54) 

T 

One method to allow for relaxation of the u-velocity component is to 
introduce the e_ equation, the last row of Equation (13) . Further, at 
the conical stagnation points e_ is related to p_ by 

e_ = p_ (at (|) - 0, 7r) (55) 

T YP T 

With e known, u follows from Equation (3) . 

Inflow and Outflow Meridional Boundaries 

The three subproblems making up the partitioned half-plane (Figure 2a) 
are joined by meridional inflow and outflow boundaries. 


In the windward region the outflow boundary must be selected a priori 
such that it lies beyond the limiting characteristic in the crossflow plane 
Extrapolation then provides flow values along the boundary. Where possible 
second-order, equally-^-spaced extrapolation is used. However, since the 
floating bow shock may cut across the reference mesh it is necessary to 
also allow for unequally spaced data obtained by interpolating bow shock 
values. 

The shoulder region (f-marching code (Appendix A) obtains initial data 
from the converged windward solution. The w~velocity must be supersonic 
all along the initial data line. The i|)-direction sweep around the cone is 
terminated when the w- velocity drops to sonic speed. In practice, a cut- 
off of w = 1.05a is used since the (|)-step size approaches zero as w 
approaches a. 

The leeward region inflow boundary flow variables are held fixed at 
values prescribed by the shoulder region code. The initial meridional 
plane is selected such that the crossflow Mach number always exceeds unity. 

Symmetry Planes 

Imposing the symmetry conditions. Equations (53), on the first four 
rows of Equation (13) yields the three equations to be integrated at the 
windward and leeward symmetry planes. The symmetry conditions are also 
applied in the calculation of the vortical singularity and in the calcula- 
tion of the body and bow shock points at the symmetry planes. 

As with the w- derivative in Equation (52) , the differencing of v- 
Y X 

in the v_ equation should be made in accordance with the sign of the 
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v-velocity (19) . Additionally, special discretization is required at 
points neighboring the vortical singularity in the leeward symmetry plane. 

The updated entropy is specified above and below the vortical singu- 
larity and e computed as in the body boundary calculation. 

Symmetry conditions may also be applied to the conservation-law form 
of the governing equations (Equation (9) ) . If the conservation- law form is 
integrated in the leeward symmetry plane, no special treatment is given to 
the vortical singularity. 

Discontinuities 

The bow shock in the shock-as-a-boundary code is calculated following 
Thomas' "pressure approach" (48) . Its adaptation to the shoulder region 
problem is detailed in Appendix A. In the windward, leeward, and symmetri- 
cal half-plane problems the use of Thomas’ scheme is similar to the appli- 
cation discussed in Reference 31, 

The floating-fitting technique utilizes Kentzer's characteristics 
based approach (52) to propagate the bow shock, the embedded crossflow 
shock, and the vortical singularity. 

Detection and monitoring 

Embedded discontinuities are searched for in the flowfield provided 
by the initial condition procedure (Chapter IV). The floating vortical 
singularity is initially positioned at a weighted location within the mesh 
interval having the maximum density gradient in the leeward symmetry plane. 
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Crossflow shock points are detected by continually scanning the 
pressure distribution along 6 = constant lines. Trial floating shock 
points are positioned at a weighted location within the mesh interval 
having the maximum compressive gradient (as computed by the conservation 
law form of the governing equations) . With extrapolated upstream flow 
values and a finite-difference approximation for the shock slope, a 
normal Mach number is computed. If the normal Mach number is greater 
than one, then tracking of the trail shock point as a shock point 
begins. 

A check for embedded shock points is made after each advancement, 
with points being added or discarded. The embedded shock forms early 
and locks into place near the cone surface. However, the end-point of 
the embedded shock tends to oscillate. In order to eliminate this 
problem, without the complications of treating the true shock tip 
(41), an artificial cutoff of the fitted shock is made. Some small 
overshoots and undershoots, ti'pical of shock-captured solutions, are 
observed at the truncated shock tip. However, these small errors do 
not propagate away from the tip region. 

In practice , the simplest of several schemes used to terminate the 
crossflow shock is to specify a priori a 9-boundary beyond which the 
crossflow shock is captured. 
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Propagation 

Combining a characteristic compatibility relation with the differen- 
tiated discontinuity j\amp conditions yields an equation for the accelera- 
tion of the discontinuity. The updated discontinuity speed and position is 
obtained by integrating its acceleration using a second-order Euler pre- 
dictor/modified Euler corrector method. 

The updated flow values on one side of the discontinuity are obtained 
by integrating the non conservative form of the discontinuity aligned gov- 
erning equations (Equation (13) ) , with the exception of the bow shock where 
one side is the known free stream. At the crossflow shock, the governing 
equations are integrated on the low pressure (crossflow upstream) side. 

At the vortical singularity, the integration is performed on the low en- 
tropy (bow shock) side. 

With the updated geometry and flow values known on one side , the up- 
dated flow values on the opposite side of the discontinuity are determined 
by the jump conditions. 

The implementation of Kentzer's scheme may be illustrated by consider- 
ing the propagation of the crossflow shock. Crossflow shock points are 
tracked in the Y-direction. Using the notation and definitions shown in 

Figure 6 , q^_ represents the acceleration of the crossflow shock, normal 

■ . . 

to the shock and pointing towards the high pressure side. The equation for 
the crossflow shock acceleration (with clustering in the reference mesh) 
obtained by combining the compatibility relation along the down-running 
characteristic in the Y 3 , T-plane (Table 1) with the f -differentiated 
Rankine-Hugoniot jump conditions, is 
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In Equation (56) , the down-running characteristic is defined by 
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the upstream normal relative Mach number is 

r>j 
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and the nonsubscripted flow values represent the crossflow shock downstream 
conditions . 

The crossflow shock acceleration (Equation (56)) is used to update the 
normal component of the shock velocity and the shock position according to 






the two-step sequence 


Predictor: 


n+1 n . n . _ 

q~ = q^ + q^ At 

s ^s ^s_ 

T 

,n+l j.n ,n 

9~ = At 

s ^s s_ 

T 


Corrector: 


n+1 n . 1 n+1 . n 

% 2 

V T T-* 

,n+l n . 1 f.n+1 , .n 1 
s 2 ^s_ s_ 

V. T T' 


As is shown in Figure 6 , the component of the shock velocity in the 
9“direction is related to the normal component by 




‘f’^T /^sq'^ sin 9 *^3 


With the updated crossflow shock low pressure side flow values, the 
upstream velocity component normal to the crossflow shock is given by 


ui = VI ng + wi n^ 


and, relative to a shock moving with velocity q^. 
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Application of the Rankine-Hugoniot relations yields 
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The u-velocity component (velocity in R-direction) is tangential to 
the crossflow shock and hence remains unchanged 

U2 = ui (69) 


The remaining downstream velocity components are 


V2 = Ilg + VI 


W2 - + W1 


(70) 

(71) 


with the velocity components determined, the total energy is obtained from 
Equation (3) . 

The propagation of the bow shock is analogous to, but much simpler 
than, that described for the embedded crossflow shock. 

The motion of the vortical singularity in the leeward symmetry plan^ 
is, following Kentzer's approach, governed by 




{[(ai +a 2 )p^ + YP(v 2 y- vi-)] x 
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where is the speed of the vortical singularity, positive pointing 

towards the body, and where the subscripts 1 and 2 represent, respec- 
tively, the bow shock and body sides of the singularity. 

Prior to advancing the vortical singularity, the bow shock in the lee- 
ward symmetry plane and the crossflow shock at the cone surface are advanced 
to make available the updated values of leeward symmetry plane entropy. 


denoted as Si and S 2 . 
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With the updated position and speed of the vortical singularity, to- 
gether with the bow shock side (side 1 ) flow values, the body side pressure 
is given by rhe continuity of pressure across the singularity 

P2 = Pi (73) 

and the absence of a normal relative velocity component yields the body 
side v-velocity 

V2 = VI = - (74) 

Equating total enthalpy across the singularity yields an equation for 
the R-direction velocity component 
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then Equation (3) is used to compute eg • 

Alternatively, the total enthalpy need not be equated across the sing- 
ularity provided the e- or u- equation is integrated on the body side of 
the singularity. 

Variable area effect 

As noted in References 39 and 41, shock acceleration equations can be 
very sensitive to the manner in which they are discretized. The crossflow 
shock acceleration Equation (56) is no exception. 

Following the suggestions in Reference 39, the R terms in Equation 

(56) have been grouped according to their differing physical roles. The 

accurate calculation of R 5 requires that the v- and w- difference ap- 

X X 

proximations include discontinuity mesh points which lie midway between the 
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mesh points being tracked. The inclusion of this so-called "variable-area" 
effect thus doubles the niomber of discontinuity mesh points. The location 
of these intermediate mesh points is determined by the points being tracked. 
The evaluation of side 1 temporal derivatives at the intermediate mesh 
points introduces the further complication of obtaining upstream flow values 
by interpolation in the reference mesh. However, with these added mesh 
points the tendency of the crossflow shock to develop kinks is removed. 
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CHAPTER VI. RESULTS AND DISCUSSION 

The partitioning of the crossflow plane into separate windward, shoul- 
der, and leeward region problems (Figure 2a) is demonstrated by solving the 
case: = 7, 9^^ 20°, and a = 30°. With this choice of parameters, re- 

sults of the floating-fitting method can be compared with those obtained in 
Reference 34 using the GTT modified method of lines approach. 

Additional windward region results are obtained for the case : M - 

OO 

7.95, 0^ = 10°, and a = 16° to compare the numerical floating-fitting solu- 
tion with experimental measurements reported in Reference 47. The compari- 
son is made on the windward side of the cone where viscous effects are 
small. Experimental data on the leeward side cannot be used to verify the 
numerical method since discrepancies would be due primarily to the inappro- 
priateness of the inviscid model to describe the viscous dominated region 
near the cone surface and in the symmetry plane. 

The symmetrical half -plane problem (Figure 2b) , where partitioning of 

the crossflow plane is not possible, is demonstrated by computing the case: 

M =3, 6 - 1 . 5 °, and a = 15°. In Reference 31, a solution obtained using 

00 c 

a finite-difference method based on the shock-capturing approach (for 
embedded discontinuities) is available for comparison. 

The stability and accuracy of the floating-fitting technique was first 
tested on the windward region problem, thus avoiding the additional com- 
plexity of embedded discontinuities. Also, several boundary condition 
schemes were assessed using windward region calculations. 
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Figure 7 shows the solution for the floated bow shock shape in the 
reference mesh compared with the shock-as-a-boundary r’esults and a method 
of lines solution (34) . The corresponding density distribution behind the 
bow shock is plotted in Figure 8 with the method of lines solution from 
Reference 13. The floating-fitting, shock-as-a-boundary, and method of 
lines results are in excellent agreement. Numerous crossings of the refer- 
ence mesh 6 = constant grid lines (Figure 7) are demonstrated to have no 
adverse effect on the floating-fitting solution. 

Body boundary condition methods are compared in Figure 9 for a wind- 
ward region case where experimental values are also available (47). The 
numerical methods all tend to produce the same result- The differences be- 
tween the numerical and experimental pressures are primarily due to experi- 
mental errors (12, 52) with viscous effects being negligible (53). 

It might be expected that Kentzer's body boundary scheme, based on the 
theory of characteristics, will yield the most accurate solution. However, 
the windward region results in Figure 9 do not distinguish one body boun- 
dary procedure as being superior. This close agreement of the various 
techniques is most likely due to the use of coordinates and velocity com- 
ponents (independent and dependent variables) which are in the surface 
normal and tangential directions. 

In Figure 10, the solution from the shoulder region (j) -marching code 
is shown together with the windward region results in the form of crossflow 
plane contour plots. The acceleration around the cone shoulder is demon- 
strated by the crossflow Mach number levels in Figure 10a. The windward 
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crossflow sonic lines lies well upstream of the windward outflow boundary 
located at = 90° . 

The constant entropy ;’ontours in Figure lOh represent the streamline 
pattern. The shoulder regi.on contours match precisely those obtained by 
the method of characteristics (13, 34). 

Results for the leeward region flowfield are presented in Figure 11. 
The floating bow shock dips slightly in the leeward symmetry plane (Figure 
11a) . The shock shape thus represents the so-called "anomalous" position 
(23, 35) and differs from the "regialar" position of the extrapolated bow 
shock. 

The analysis presented in Reference 35 lends support to the floating- 
fitting solution for the bow shock shape. In Reference 35 an expression is 
derived, based on the assumption of small incidence, for the angle of at- 
tack at which the bow shock transitions from the regular to the anomalous 
position. For M =7,9= 20°, and y = 1.4 the maximum shock standoff 

TO Q 

distance moves out of the leeward symmetry plane when a > 4.1°. The for- 
mation of the crossflow shock and the lift off of the vortical singularity 
at large angles of attack is not accounted for in this analysis. Results 
presented in Reference 54, however, do Verify the change to the anomalous 
position for a = 5° and the movement of increasingly away from 

the symmetry plane for angles of attack up to 15° . 

The density distribution behind the bow shock is shown in Figure lib. 
The solution verifies the conjecture (13) that for relatively "thick" cones 
with "moderate" or large free-stream Mach numbers it might be expected that 
the bow shock will retain nonzero intensity in the leeward symmetry plane. 
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That is, although the bow shock dips in the leeward synunetry plane, it does 
not reach the free-stream Mach cone (Figure 11a) . The density behind the 
bow shock in Figure lib drops to about 1.2 times the free-stream value. 

The crossflow shock, in Figure 11a, is perpendicular to the surface at 
its base. Away from the surface, the crossflow shock bends very slightly 
towards -the crossflow upstream direction. The normal intersection of the 
crossflow shock, with the cone surface is not imposed in the floating- 
fitting code, but evolves as the steady state conical solution is approached 

The length of the crossflow shock in Figure 11a is the length of the 
fitted portion only. As described in Chapter V, the fitting scheme is not 
applied all the way to the crossflow shock tip (where the upstream normal 
Mach number is one) . 

The pressure and density distributions on the downstream side of the 
crossflow shock are shown in P.lgure 11c. Tentative evidence for' the pre- 
sence of a conical logarithmic singularity, analogous to that found in two- 
dimensional flow, is provided by noting the difference between the Rankine- 
Hugoniot determined p slope at the cone surface with that from the 

D ’ 

normal momentum equation with v = 0 ;; 

p. = pw^ cot 0 (at 0 = 0 ) (76) 

■0 G 

(plotted as a dashed line in Figure 11c). 

This numerical evidence is, of course, only speculative. In lieu of 
an analytical proof, further numerical parametric studies are needed to 
support this conjecture on the existence of a logarithmic singularity. 
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The leeward symmetry plane density distribution, obtained by solving 
the conservation-law form of the governing equations without any special 
treatment of the Vortical singularity, is compared with a modified GTT 
method solution (34) in Figure lid. The strong density discontinuity at 
the vortical singularity is spread out over several 0 -degrees in the 
singularity-captured solution. From the leeward symmetry plane density 
(and u-velocity) distributions it may be inferred that the singularity- 
captured solution demonstrates, at least qualitatively, lift off of the 
vortical singularity. However, the smearing out of the singularity, which 
closely resembles the viscous behavior observed experimentally, gives rise 
to the largest total enthalpy errors (‘-“0(1%)) in the flowfield (see Accu- 
racy and Convergence Section of Chapter III) . The 0-direction velocity 
component in the leeward symmetry plane (v-velocity) is negative at the bow 
shock, stagnates at the vortical singularity and again at the cone surface, 
and should be positive between the singularity and the surface , With a 
singularity-captured solution small oscillations in the v-velocity distri- 
bution between the singularity and the surface (where the magnitude of v 
is small) can cause negative v-velocities . Thus, some flov; features near 
the singularity are not resolved using the singularity-capturing approach. 

Implementation of the floating-fitting procedure for the vortical 
singularity, presented in Chapter V, does not improve the resolution of the 
singularity. The fitting technique requires accurate information from the 
flowfield interior (transmitted by the w^ derivative) which is computed 
using the conservation-law dependent variables. However, large gradients 
in the flowfield, but not discontinuities, exist adjacent to the leeward 
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symmetry plane. Thus, unlike the flow regions on each side of floating- 
fitted shocks, the neighborhood of the singularity is not smooth, especially 
in the 6 -direction. The conservation-law formulation tends to smooth these 
gradients. The floating-fitting scheme, in turn, moves the singularity and 
imposes jump conditions based on the smoothed interior flowfield values. 

The result of this complicated interplay between the singularity, the sym- 
metry plane boundary conditions, and the flowfield interior is that a small 
vortical singularity velocity remains (nonconvergence) while the neighboring 
flowfield still resembles the singularity-captured solution. A new set of 
dependent variables for the flowfield interior, similar to those suggested 
in Reference 55, might help to resolve the singularity. However, the 
conservation-law form has proven to be very convenient in regards to the 
crossflow shock end-point treatment; such benefits would be lost with a 
change of dependent variables. 

A S'ymmetrical half-plane solution is presented in Figure 12 . The bow 
shock loses its strength on the leeward side as it approaches tangency with 
the free-stream Mach cone (Figure 12a). The density behind the bow shock, 
plotted in Figure 12b, thus becomes equal to that of the free stream on the 
leeward side. Kentzer's shock boundary scheme accurately maintains the 
zero strength shock solution. 

The floating-fitting result is in agreement with the conclusion in 
Reference 16 that for relatively "thin" cones or small free-stream Mach 
numbers the bow shock will degenerate into a Mach wave in the leeward 
symmetry plane. 
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As shown in Figure 12a, the crossflow shock (fitted portion) does not 
extend very far into the flowfield interior. The pressure jump across the 
crossflow shock is apparent in Figure 12c where the surface pressure dis- 
tribution is shown. The floating-fitting results in Figure 12c match those 
obtained by a shock-capturing approach with slight differences occurring 
downstream of the crossflow shock. The captured internal shock appears as 
a sharp jump due to the alignment of the shock with the mesh, the clustering 
of points near shock, and its relatively weak strength. 
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CHAPTER VII. CONCLUDING REMARKS 

Complicated flowfields occur about the geometrically simple circular 
cone at large angles of attack in a supersonic free stream. The cone prob- 
lem serves as convenient test of n\americal techniques intended for the com- 
putation of multidimensional flowfields containing embedded discontinuities. 
The conical nature of the problem provides a clear accuracy assessment that 
is not available in the calculation of flows about complex nonconical con- 
figurations. 

All of the established numerical methods have encountered difficulties 
with cone at large angle of attack calculations. The most successful tech- 
nique to date, a modification of the method of lines, is restricted to cases 
in which the crossflow sonic lines extend to the bow shock. 

In this study an explicit finite-difference technique, based on the 
concept of floating- fitting discon inuities, has been described. Results 
have been presented which demonstrate the capabi3.ity of the method over a 
range of Mach numbers, cone angles, and angles of attack. The technique 
has been shown to accurately compute both strong embedded shocks and van- 
ishingly weak peripheral shocks. The method does not resolve the structure 
of the flowfield in the immediate vicinity of the vortical singularity. 

From a pragmatic standpoint, the computer code for the floating-fitting 
algorithm, which must include logic for numerous discontinuity crossings of 
the computational mesh, is necessarily lengthy. However, the execution 
tiraes are not excessive since most mesh points require no special treatment . 
The necessity to input several a priori estimates and generate initial con- 
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ditions coupled with the overall complexity of the code makes each case a 
trial and error process. 

From a theoretical viewpoint, it is recommended that future efforts 
might be directed towards refining the shock tip calculation. Also, the 
treatment of the vortical singularity might benefit from a study of the 
mutual Influence of boundary condition schemes with various sets of depen- 
dent variables in the flowfield interior. 
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APPENDIX A. SHOULDER REGION CALCULATION 


The shoulder region flowfield about highly inclined cones is most ef- 
ficiently calculated by a two-dimensional meridional marching algorithm 
rather than a three-dimensional iterative relaxation approach. The proce- 
dure requires that the parameters a, and 0^ are such that there exists 

a region- of supersonic meridional velocity component extending from the cone 

surface to the bow shock (Figure 2a) . The principle features of the tech- 
's 

nique are outlined below. 


Dependent and Independent Variables 

The steady Euler equations in cylindrical polar coordinates (z,r,<j>) 
with the z-axis coincident with the cone centerline, r the cylindrical 
radial direction, and <f) the meridional angle measured from the windward 
symmetry plane, may be written in weak conservation-lav/ form as 

E' + F’ + G! + H‘ = 0 (Al) 

zr <p 

where E', F', G' , and H' are the four-component vectors 




F’ 


pv 

puv 
— • ~ 


kp + pv' 
pvw 





pw 

pwu 

pwv 

kp + pw^ 



pv 

]^,puv 

r p (v^ - w^) 
2pvw 


(A2) 


The velocity components u, v, and w are in, respectively, the z, r, 
(f)-directions . The pressure p and density p are nondimensionalized with 
respect to their free-stream stagnation values and the reference speed is 
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the maximum adiabatic speed. The constant k is related to the ratio of 
specific heats y by 


k = 


Y - 1 

2y 


(A3) 


The governing equations are made complete by the steady flow conserva- 
tion of energy relation 

p = p (1 - u^ - v^ - w^ ) (A4) 

In order to treat the bow shock as an outer boundary of the computa- 
tional domain, the independent variable transformation, 


C = z 

r - r (z) 

. ^ 2 

r^((|), z) - r^(z) 

$ = (f) 


(A5) 


v;here 

r = z tan 0 ~ cone radius 

c c 

r {(|), z) ~ bow shock radius 
s 

is applied to Equation (Al) . With the assumption of conical flow. Equation 
(A5) represents a self-similar transformation and all flow derivatives with 
respect to C vanish. The transformed conical governing equations are 

G| + F* + H* = 0 (A6) 


where 
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with 


G* = G' 

F* = E' + F' + G' 

H* = H' - E’ - P' - G’ {^)^ 



(A7) 


(A8) 


Finite-Difference Method 

With initial conditions specified along a (}> = constant line (provided 
by a windward region solution), MacCormack's predictor-corrector algorithm 
(Appendix B) is applied to Equation (A6) to advance the solution in the 
meridional direction. The calculation is terminated when the minimum Mach 
number in the hyperbolic marching direction falls below 1.05. 

The F* derivative at interior mesh points is approximated using the 
forward-predictor, backward-corrector version of MacCormack ' s algorithm. 

At the cone surface , F| is always approximated by forward differences 
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(two-point) with all dependent variaCles being overwritten by the body 
boundary condition procedure following the corrector step. At the bow shock 
boundary, F* is calculated using backward differences (two-point) in both 
the predictor and corrector steps. With given free-stream conditions, in 
conjunction with a sharp shock boundary procedure, information is trans- 
mitted across the bow shock only through the Rankine-Hugoniot relations. 

The shock boundary condition scheme is applied following both the predictor 
and corrector steps. The finite-difference value for pressure is accepted 
with the remaining flow variables being overwritten by the sharp shock pro- 
cedure . 

Stability Analysis 

An estimate of a stable step size for the explicit MacCormack algorithm 
is provided by the C-F-L (Courant-Friedrichs-Lewy) condition 




(CN)Ag 


max 


(A9) 


where 


r/ i w (au f v) + (w~ - a^) (b/r) j + a/ (au + v) ^ + (w^ - a^) (a^ + 1) 


max ' 


,,2 i .2 

W — a. 


(AlO) 


represents the slope of a characteristic projected onto the 5, #-plane 
with 


b = 


(All) 


c 
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and with the local speed of sound 



(A12) 


From linear theory, the Courant parameter CN should be less than or equal 
to 1.0; typically, for shoulder region calculations CN is set to 0.9. 


Surface Boundary 

Surface tangency, constant total enthalpy, and known surface entropy 
are imposed according to Abbett's (50) body boundary scheme. 

Let a superscript (1) denote the surface flow values obtained from the 
finite-difference algorithm following the corrector step. In order for the 
velocity vector 


=u'i>£ 

z r (j) 


(A13) 


to satisfy surface tangency, the flow is turned through a small isentropic 
compression or expansion where the turning angle is 


Av = sin“^ 


C J- ) A 

q • n. 


( 1 ) 


(A14) 


with 


n 


1 +1 
^Z Z 2 


(A15) 


the unit vector normal to the cone surface. 

The series relations for the pressure following a turn through an 
oblique shock wave or an expansion through the small angle Av are 
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equivalent if terms O(Av^) are neglected. With Av from Equation (Al4) 
and the Mach number 


M 


M 



( 1 ) 

IT) 


(A16) 


the pressure resulting from the tangency requirement is approximated by 


p ^ 1 - 


Av (y + Dm”^ - 4 (M^ - 1) 


/m2 - 1 


4 (m2 - 1) ^ 


(A17) 


Since surface entropy is specified in the initial meridional plane and 
remains constant along the body streamline, the surface density follows 
from 




(A18) 


The velocity modulus satisfying constant total enthalpy is obtained 
through the energy relation (Equation (A4)) 


q = t^l - p/p' 


(A19) 


With the magnitude q and direction t^ known (where t^ denotes a unit 
vector tangent to the cone surface) , the finite-difference velocity com- 
ponents u^^\ v^^^ , and w^^^ are replaced by 


u = — ^ 


[tI 


u''^^ + N 




q (1) 
w = ' ^ - w 


(A20) 



where N is the scalar 


^(1) 
q • n 


N = 


Irl +1 

'-z 


( 1 ) ^ ( 1 ) 

-u r„ + V 


(A21 


+ 1 


and T the tangent vector 


T = q 

^( 1 ) 


— ( 1 ) <> (A Irpia- 


- N 


-TZn 1+1 


'z z r 


(A22: 


with magnitude 


T = 


[ (1) , „ 1 

2 

f (1) ' 

2 

\ (1)^ 


u + N r- 

1 

+ 

V - N 

+ 

w 



(A23; 


Shock Boundary 


The two-dimensional shoulder region calculation is performed in the 
plane perpendicular to the cylindrical z-axis. Since the flowfield is 
self-similar, the slope of the bow shock in the z-direction is simply 


s 

z 


(A24: 


The slope of the bow shock in the marching meridional direction rg^ 
is obtained by inverting the relation for the component of the free-stream 
velocity normal to the bow shock ui , where the subscript 1 refers to 
the bow shock upstream conditions. The resulting two roots are 
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Wj < ± Uj 



(A25) 


where 


^3, 


K = 


2 ~2 
- "I 


(A26) 


with the free-stream velocity components 

ui = qi cos a 
■'^1 “ “1l sin « cos (|) *■ (A27) 

wi = sin a sin cf) ^ 

In Equation (A25) the appropriate root when w| > contains the + sign 
while the - sign is used if < u^. The exact equality Wj = u^, v/hich 

would lead to a zero in the denominator of k, is not encountered numeri- 
cally. In typical shoulder region calculations the switch in sign occurs 
without any difficulties at those meridional stations bounding the (j)- 
location at which Wj == u^. 


The bow shock computational plane boundary in the shoulder region code 
is computed following a "pressure approach" (48) . With the known initial 
bow shock downstream pressure P 2 (where P 2 = P2 prior to the predictor 


n+1 


step and P 2 = P2 to the corrector step) , the corresponding up- 

stream normal velocity is given by the shock jump condition 




V - 1 PI P2 I If - 1 
Pi [pi Y + ij 


(A28) 
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Substitution of Uj into Equation (A25) gives the rate of change of 
the bow shock radius in the inarching direction. The bow shock is advanced 
using the Euler predictor 


n+1 

r 

s 


n 

r 

s 


+ A(j) rc 


(A29) 


and modified Euler corrector 


n . n+1 


sequence . 



(A30) 


At the advanced (})-level the pressure behind the bow shock, obtained 
from the MacCormack algorithm with one-sided differencing in the ^-direc- 
tion, is the only flow variable accepted. The finite difference density 
result is discarded in favor of the Rankine-Hugoniot value 



P2 




Y -1 

Y + 1 


1 + '^“^ 

P2 

Y + 1 



(A31) 


Equation (A28) gives the updated Uj and Equation (A25) supplies the up- 
dated shock slope “ Finally, the updated velocity components, computed 

from the shock jump relations applied at the updated shock position, are 
(the shock having been propagated by Equation (A29) or Equation (A30) ) 


ui I (1 -P1/P2) 

V£ = v^ + ^ — 

U2 = U| - rg^(v2 - vi) 

^S(j) 

= wi - (V2 - Vj) 

s 


(A32) 

(A3 3) 
(A3 4) 
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APPENDIX B. INTEGRATION ALGORITHM 


The gas dynamic equations written in nonconservation-law form as 


+ tBi]d + [B2]d + A3 = 0 


(Bl) 


or cast in weak conservation- law form as 


U + F, + G, + H = 0 
T X Y 


(B2) 


are integrated at each point of the finite-difference mesh by the explicit^ 
two-step Euler predictor/modified Euler corrector sequence 


Predictor : 


^n+1 _n . , ^n 
f, „ = f + At f 
J,K T 


(B3) 


Corrector: 


pn+1 


^n+1 , .cU , . _n+l 
f + f + At f 

T 


(B4) 


where f represents either the d or U vector. The mesh point indices 
(J,K) for, respectively, the X- and Y-directions have been omitted from 


the right side of Equations (B3) and (B4) . The superscripts n, n+1, and 

n+1 refer to in order the current, predicted, and advanced time levels 
, n+1 n , . ' 

The MacCormack (49) variant of the Lax-Wendroff method, in which spa- 
tial derivatives are approximated by equally spaced, two-point, one-sided 
differences of altering direction in the predictor and corrector stages, is 
used at all mesh points not lying on boundaries or neighboring (in time or 
space) discontinuities. The resulting noncentered algorithm is second- 
order in both time and space . 


For all floating-fitting applications forward differencing in the pre- 
dictor is followed by backward differencing in the corrector (several per- 
mutations of the differencing direction sequence were tested and found to 
yield substantially the same results for the windward region of the flow- 
field about a cone at large angle of attack) . 


To illustrate; the MacCormack algorithm applied to Equation (B2) to 


- H 


advance the 

point (J,K) 

to time leve 

n+1 

T 

approximates 

(B3) 

by 



F - F 

J+1 J 


r' _ r* 




_ 

AX 

AY 

W 4 . kjl 

_ K+l K 


and 

T 

in Equation 

(B4) by 







_ J:- 

T? _ 1? 

_ A- 

r’ ... r* 




AX 


AY 

\ ViJ 



T 


(B5) 


- H 


(B6) 


where the nonvarying sub- and superscripts on the right sides of Equations 
(B5) and (B6) have been deleted. 

On boundary points , where only one direction is available for forming 
differences in both the predictor and corrector stages, equally-spaced 
three-point approximations (which are consistent with the overall second- 
order of the scheme) are used. 

In the presence of discontinuities the usual MacCormack algorithm is 
modified to prevent forming differences in time or space that would cross a 
discontinuity surface. Where required, unequally spaced one-sided spatial 
derivative approximations replace the approximations in Equations (B5) and 
(B6) . These modifications to MacCormack ' s scheme, necessitated by a 
floating-fitting approach, are dealt with in Chapter III. 
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